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A simple procedure to approximate the noncommutation terms that arise whenever it is necessary 
to use a variable scale filtering of the motion equations and to compensate directly the flow solutions 
from the commutation error is here presented. Such a situation usually concerns large eddy simu- 
lation of nonhomogeneous turbulent flows. The noncommutation of the average and differentiation 
operations leads to nonhomogeneous terms in the motion equations, that act as source terms of in- 
tensity which depend on the gradient of the filter scale 8 and which, if neglected, induce a systematic 
error throughout the solution. Here the different noncommutation terms of the motion equation 
are determined as functions of the 8 gradient and of the 8 derivatives of the filtered variables. It is 
shown here that approximated noncommutation terms of the fourth order of accuracy, with respect 
to the filtering scale, can be obtained using series expansions in the filter width of approximations 
based on finite differences and introducing successive levels of filtering, which makes it suitable to 
use in conjunction with dynamic or mixed subgrid models. The procedure operates in a way which 
is independent of the type of filter in use and without increasing the differential order of the equa- 
tions, which, on the contrary, would require additional boundary conditions. It is not necessary to 
introduce a mapping function of the nonuniform grid in the physical domain into a uniform grid 
in an infinite domain. A priori tests on the turbulent channel flow (Re T = 180 and 590) highlight 
the approximation capability of the present procedure. A numerical example is given, which draws 
attention to the nonlocal effects on the solution due to the lack of noncommutation terms in the 
motion equation and to the efficiency of the present procedure in reducing the commutation error 
on the solution. 



I. INTRODUCTION 

The problem of the non commutativity of the fil- 
tering operation has been considered by Ghosal and 
Moin (1995)Q, van der Ven (1995)0, Fureby and Ta- 
bor (1997)[|, Vasilyev et al. (1998) [|] and reviewed by 
Ghosal (1999)@ and Sagaut in his monography on large 
eddy simulation (LES) for incompressible flows (2001) [6]. 

Ghosal and Moin[l[ showed that the commutation er- 
ror is of the second order in the filter width. Introducing 
a filter definition built on the mapping function of the non 
uniform grid, they proposed a procedure that can be used 
whenever numerical schemes based on pseudo-spectral 
methods or on finite differencing of an order higher than 
the second order are employed. The non commutation 
term are expanded in a Taylor series in the filter width, 
where the coefficients depend on the spatial derivatives 
of the filtered field and the mapping function. In this 
way, extra terms appear in the filtered equations, which 
increase the differential order of the equations. The au- 
thors suggested both the use of additional boundary con- 
ditions to mantain the well-posedeness of the problem or 
the use of asymptotic expansions of the filtered variables, 
in terms of the square of the filter width, which however 
requires the solution of additional non homogeneous per- 
turbative problems. 

A family of one parameter filters commuting with dif- 



ferentiation up to any given order in the filter width, 
which is assumed nonuniform in the integration domain, 
was constructed by van der Ven (2j. If a discretization 
scheme of a given order is adopted in a LES, one may 
select a filter inside the family so that the lack of com- 
mutation between differentiation and filtering can be ne- 
glected. 

A general formulation of the commutation error, due to 
a non uniform filter and to the presence of boundaries, 
has been proposed by Fureby and Tabor Boundary 
domain terms are explicitly formulated inside this repre- 
sentation. A detailed numerical analysis of the field dis- 
tribution of the intensity of the non commutation terms is 
given in this paper by comparing LES and direct numer- 
ical simulation (DNS) data obtained from simulations of 
the incompressible turbulent channel flow at Re T — 180 
and 395. It has been found that the local intensity can 
be as high as 21% of the local advection, with a field 
volume averaged relative intensity of about 8%. Lower 
values apply if reference is made to the sum of the local 
advection, pressure gradient, and molecular viscous flux 
terms. An interesting result is that the use of different 
subgrid scale models negligibly affects the local and av- 
erage values of the sum of the commutation error terms. 
As expected, the high relative intensities of the commu- 
tation terms are concentrated in the flow regions where 
the gradient of 8 is high. However, the question of the 
possible extension of the effects on the macroscopic scale 
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of the flow is still left open. 

A generalization of the procedures proposed by Ghosal 
and Moin and van der Ven is presented in the paper by 
Vasilyev et aZ.Q. A minimization of the commutation 
error is achieved by using a class of filters with n — 1 van- 
ishing moments, where n is the order of the employed nu- 
merical discretization scheme. The authors also supply a 
group of rules to construct discrete filters that commutes 
with differentiation up to any given order inside complex 
domains. 

The method here proposed relies on an approximation 
of the specific non commutation term that corresponds to 
the different terms of the motion equations. A commuta- 
tion approximation of the fourth order in the filter width 
can be obtained thanks to the introduction of successive 
levels of average. See Sec. |TT] for the basic formulation 
relevant to an isotropic grid stretching and the Appendix 
for the more general anisotropic case, which also speci- 
fies the formulation that is appropriate to wall-bounded 
flows. 

While performing large eddy simulations, the present 
approach can conveniently be used together with sub- 
grid models based on analog multi-level filtering, e.g. 
models which apply the dynamic procedure, Germano 
et al. (1991)0, Germano (1992)@ or the Bardina mixed 
model (1980) [9]. The filtering approach we use in this 
paper is that of the very fundamental volume average, 
first applied by Smagorinsky in 1963 10]. The volume 
average formulation is advantageous because it does not 
introduce an error associated to domain boundaries, thus 
avoiding the problem of the addition of further non com- 
mutation terms in the equations. 

The variable scale filtered Navier-Stokes equations, in- 
cluding the commutation terms approximation, are given 
in Sec. flTAl 

A priori tests on the turbulent channel flow data bases 
by Alfonsi et al. (1998) [llj and Passoni et al. (1999) [HI, 
Re T = 180, and by Moser et al. (1999)0. Re T = 590, 
are presented in Sec. IIII Al An example of application of 
the numerical procedure is presented in Sec. IIII B[ which 
focuses attention on the fact that this systematic error 
is important throughout the entire flow and not only in 
the regions where the non homogeneous terms of the mo- 
tion equation, which originate from the lack of commuta- 
tion of the operations of differentiation and filtering, are 
different from zero. The capacity of the present proce- 
dure to reduce the relevant absolute and relative errors 
is shown. 

Before proceeding to the other sections, it is necessary 
to open a digression on the terminology adopted in what 
follows. Since, among the points being discussed in the 
paper, there are: the structure of the motion equation, 
once a variable scale filtering is used, and the role played 
by the terms which originate from the non commutation 
filter-differentiation, we have had to clearly distinguish 
the concept of commutation error on the flow solution 
from that of non commutation term in the equations. 
By commutation error we mean the error which affects 



the flow solution, when a variable scale filter is used, but 
the equations are used as if the operations of filtering and 
differentiation commute. By non commutation term we 
mean any of the terms which originate in the equation 
of motion when the filter scale is a function of the point. 
It is necessary to recall that in previous literature the 
latter was called commutation error, since the equation 
were always used as if they were commutative and the 
omissiom of the terms, which should make them complete 
when the filter length varies, introduced the error in the 
solution. 



II. NON COMMUTATION TERMS AND THEIR 
APPROXIMATION 

The loss of the commutation between the spatial fil- 
tering and the differentiation operations is related to the 
use of a variable filter, which in the more general config- 
uration is anisotropic. In this case the filter width is the 
vector <5(x) = (5i(x), S 2 (x), <5 3 (x)). 

For reading convenience and as the isotropic stretching 
configuration is conceptually non reductive, a scalar filter 
scale <5(x) is assumed in what follows. However, the gen- 
eral anisotropic configuration of stretching is dealt with 
in the Appendix, which also specifies a filtering formula- 
tion that is suitable for wall-bounded flows. 

Let us suppose we have chosen a given class of integra- 
tion volumes 

V S = {r]elR 3 :|| t] || (6} 
and an average operation for the variable /(x): 

(fU = ^[ f(x + v)dv=^[ /(x + 5€)d€, (1) 
vs Jv s vi Jv 1 

where the transformation £ = rj/S has been used and, as 
a consequence, V\ — V$/5 3 . Please note that, with this 
choice, the width of the averaging volumes is twice the 
filter scale. 

A variable filter scale is introduced by allowing 8 to be 
a function of point, S — S(x). In this case 



= <V/-«i)« + ^;(x)(V/-0 
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of the derivative is a differential operator acting on the 
filtered field: 
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The non commutation term C[, which is defined as 



(4) 



(5) 



can be represented through @ by the product of the filter 
space derivative and the filter derivative of the filtered 
variable: 
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The here proposed method is based on an approxima- 
tion of relation ©. The problem could be faced adopting 
a truncated series expansion of (f)g in terms of powers 
of 8 (Ghosal and Moin, 1995) [lj. However, this would 
increase the order of the equations, and thus require ad- 
ditional boundary conditions. Here a numerical approx- 
imation of the 8 first derivative is used in conjunction 
with truncated 8 expansions. Let us write the second 
order finite difference approximation 



d(f)s _ 1 
88 2h 

Choosing h = 8 



((f)s +h -(.f)s- h ) + 0(h 2 ) 



(7) 



^ i = ^((fhs-(f}o) + 0(S 2 ) (8) 

Now, the problem to face is that of the approximation of 
(/)o = / an d (f)2S in terms of relevant averaged quan- 
tities. Using a Taylor expansion of the integrating func- 
tion in (TTJ, we obtain the following expression for (f)s, 
in terms of /, and the filter width: 

(/) 5 (x) = /(x) + i aiAO V 2 /(x)<5 2 + 
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h,o,o(^ + 9 2 4 + 9 3 4 )/(x) 



+ 6a hl , (d 2 d 2 + 8l8l + 8 2 2 8 2 3 )f(x)] <5 4 + 0(8 e ) 
^f^)+F 1 [f]8 2 + F 2 [f}8 4 + 0(8 6 ) (9) 

where coefficients a,j/- are defined as 



and the operators Fi,F 2 as: 
^[•] = 501,0,0 V 2 . 

F 2 [-} = - [a 2 ,oAdi + d% + 81) • +6a 1 ,i, (9 1 2 9 2 2 



dldl + dldl). 



(10) 



(11) 
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From ((9]) it follows that 

f = (f)s-F 1 [f]8 2 + 0(8 i ) (13) 

and then, averaging (| 13[) on a volume of linear dimension 
28, 

(fhs = ((f)shs {Fi[f]S 2 ) 25 + 0(8*). (14) 
However, from ([9]) it can be observed that 

{Fi[f]8 2 ) 2S = Fx [f}8 2 + AFi [Fx [f}8 2 }8 2 + ... = 
= 8 2 F 1 [f] + 0(8 4 ) 

so that 

(f)28 = ((.f) S )28-F 1 [f}S 2 + 0(8 4 ) (15) 

When the expressions (| 13[) for / and (fl5| for (f)2S are 
introduced into ([8]), one obtains 

^§ l -^(((f)s)2S-(f)s) + 0(S 2 ) (16) 

When using (fT6|) . the non commutation term C[ [see ((6])] 
can be approximated by 

c' i ((M = -^^(((f}6U-{fh) (17) 



which implies 
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Cm)s)=C' i ((f)s) + g^O(8 2 ). (18) 

In order to analize the approximation error (8 i 8)0(S 2 ) 
and give a true estimate of it, let us write: 



6(x) = A<p(x), 



(19) 



where A is a reference value of the filter width which 
is usually associated to the portion of the domain where 
conditions of near homogeneity of the flow hold. Function 
tp(x) is a positive non dimensional function which belongs 
to the interval [<5 mira /A, 1], in the homogeneous region 
of the flow <£>(x) is constant and equal to 1. Function 
ip(x) varies in the inhomogeneous regions, though it keeps 
values that are greater than <5 m j n /A, which is a value that 
must correspond: (i) - when the local scale invariancc 
may be supposed - to a convenient minimum value of the 
filter width still inside the inertial range, and (ii) - when 
the local scale invariance does not hold - to a scale of the 
order of the scale which characterises the local turbulence 
structure, as, in case of wall flows, is the scale of the 
quasi-streamwise vortices peculiar to the viscous sublayer 
(see Moin and Kim, 19820, and Ghosal, 19990). 

Introducing ([TP"]) and dg = A~ 1 8 lp into and deduc- 
ing the 5 derivative from expansion (J9)) , written up to the 
fourth order of accuracy, the non commutation term can 
be estimated 

+2^ 2 (x)F 2 [/]A 4 + 0(A 6 )) (20) 
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Using expantion © twice, after having inserted (fT9|) . ap- 
proximation (|17l) can be estimated as 



-F 1 [<^ 2 F 1 [/]])A 4 + 0(A 6 )) 



(21) 



To keep the validity of estimates (f20[) and (|2lj) . care must 
be taken to select a function y(x) which, in the region of 
filter variation, also posseses first and second derivatives 
of O(l). Possible examples are trascendental functions 
such as arctan(x), tanh(x). 
A comparison of (|20[) and (f2"Tj) yields 



C' l ((f)s)-C' l ({f)s)*0(A i ). 



(22) 



Consequently, when (|17p is used, a fourth order non 
commutation term in Q is produced instead of the sec- 
ond order error, which would be obtained by totally ne- 
glecting the lack of commutation [l| . Introducing a finite 
difference approximation of a higher order than ([7]) , and, 
consequently, further levels of average, it could be possi- 
ble to increase the accuracy of the approximation of the 
non commutation term, leading to a higher order error 
in (f22|) . 

This analysis pertains to differential operators of the 
first order. The analysis is similar for the second order 
differential operators. The structure of the correction 
terms remains the same and to reach the fourth order of 
accuracy the same number of levels of average must be 
mantained. The approximation of the non commutation 
term now includes the filter of the variable spatial first 
derivatives. 

The non commutation term of the second derivatives, 
being defined by 



(23) 



can be obtained by taking the derivative of the first 
derivative ([2]) as 
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that is, 
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and can consequently be approximated using the finite 
difference for the ^-derivatives, as performed for the non 
commutation term of the first derivatives. The use of the 
standard three-point formula for the second derivative of 
{f)s with respect to 8 and of relation (fl6| for d 2 x .{f)$ 
[second term on the right hand side of equation |25p]. 
yields 



di8(x) 
' 8{x) 

{di8{x)) 2 +8{x)d 2 8{x) 
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(26) 



whose order of accuracy can be determined by again us- 
ing the trasformation 5(x) = <p(x)A. If the expansions 
of the non commutation terms C" is compared with its 
approximantion C" , it is seen that, also in this case, the 
error is 0(A 4 ). 

Even if the present analysis is based on the use of 
the volume averages, it can be observed that it remains 
valid in the case where a more general kind of filtering is 
adopted. A weight function <?(£), introduced in (p}, only 
modifies the coefficients mjk [Eq. (fTO])]. which should 
now be defined as 



j k 



However, the non commutation term (6) and its approx- 
imation remain unchanged, provided the weight function 
has a compact support. It is always possible to choose 
<5(x), so that the actual integration domain of the filter 
lies inside the flow domain (for instance setting a value 
that is lower than the distance from the wall of the first 
layer of grid points for the minimum of <5(x)). In this 
way, the compactness of the support prevents the error 
linked to the presence of finite boundaries Q. It should be 
noted that this procedure operates in the physical space 
and does not rely on the use of a mapping function of 
the non uniform grid. Centered volumes of average have 
been adopted [see Eq. JI}], even though they are not 
strictly necessary as far as the average process is con- 
sidered. However, this choice shows two advantages: (i) 
physically, when the flow is incompressible, the center of 
the volume of average is also the center of gravity and 
thus the point of application of the average momentum, 
(ii) analytically, it allows for a compact, and second or- 
der in A, representation of the non commutations terms, 
which in turn are approximated by the present procedure 
with an accuracy of the fourth order. The choice of non 
centered volume of averages, which, in principle, is math- 
ematically feasible, yields a first order in A representa- 
tion of the non commutation terms, which would also lead 
to a much more cumbersome analytical structure [l|, 0. 
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A. Non commutation terms in the averaged 
incompressible Navier-Stokes equations 

Let us consider the incompressible Navier-Stokes equa- 
tions written in the form: 

dm = (27) 
d t Ui + dj (uiUj ) + dip-v djjUi = (28) 
If a filter operator is applied the system becomes 

(d lUl ) s = (29) 

(d t Ui)s + (dj(uiUj)) s + (dip)s - v {djjU^s = 0. (30) 

When 5, the linear scale of filtering (see for instance the 
definition proposed in Sec. [TTJ) is not uniform in the flow 
domain, the averaging and differentiation operations no 
longer commute. By introducing the subgrid turbulent 
stresses R^' — (u.i)s(uj)s — (u.iUj)s and the non commu- 
tation terms C/, Ca" , for the first and second derivatives, 
as discussed in Sec. HI for the isotropic filter configura- 
tion and in the Appendix for the general anisotropic and 
the wall-bounded flow configuration [see the isotropic re- 
lations ([5]), ([6]), (|23|) and (1251) ; the anisotropic relations 
(|a3 1) - (|aT6)1 (|ATT7|, (lATTOj) , ([A~19)) ; the wall anisotropy 
relations (IA.16[) . ([A. 201) ). the averaged equations are 
written as 

(31) 



d l (u l ) 5 = -C[({u l ) s ) 



d t {ui)s + dj((ui)s{uj)g) + di(p) s -vd 2 Aui) s - djR\y 



(32) 

Together with what has been explainded in detail in Sec. 
HI] (see (fTT)) and (j2"6")l ) the present procedure approximates 
the non commutation term terms on the right hand side 
of ((311 I3"2")) with an accuracy of the fourth order. The 
correction terms can thus be represented by the follow- 
ing group of relations which are determined from the field 
information obtained through two successive average lev- 
els (the second being computed over a linear scale 28): 



djSjx) 
' 2S(x) 



[{{ui)s)2S - (ui)s] (33) 



d S(x) 

Cj(( u i)s(u 3 )s) = - 2g, x j [((Ui)d(uj)sh6 - (u t )s(u 3 )s} 

(34) 



^«P)e) = -^^[{<P)6U-<p) S ] (35) 



(d 3 S(x)) 2 + 5(x)d?6(x) 



[((Ui)s)2S " {Ui)s] 



26{x) 



R 



(37) 



(36) 



which must be accordingly modified in the case of 
anisotropy of the stretching of the computational grid, 
see the previous comments and the Appendix [ ([A.lOjl — 
dXTTb . (|AT5]) - (|AT6] l and (TCTI - dA^Ol )]. 

The adoption of the volume average allows the fil- 
tered variables to be fully supported inside the physi- 
cal domain. As a consequence, a peculiar property of 
the present procedure is that there is an absence of non 
commutation terms associated to a finite or semi-infinite 
computational domain in the filtered equations. Such 
terms arise when the filtering operator requires the ex- 
tension of the dependent variables beyond the rim of the 
domain [see Fureby and Tabor (1997) Q]- 

It can be observed that the use of the double level of 
average highlights the convenience of coupling this pro- 
cedure to subgrid models which also employ it. These 
subgrid models are the mixed model, by Bardina et al. 
(1980) 9] and, in general, all models that apply the dy- 
namical procedure, by Germano et al. (1991) [7] and Ger- 
mano (1992)0. 

The filtered equations (|2"9l I3T)]) are invariant under 
Galileian transformations. Under transformation t' — ¥ 
t, x 1 — > x + ct, a spatial variation of the filter scale S(x), 
in the x, t reference system, becomes a spatio-temporal 
variation 6(x) — 6(x' — ct') = 6(x',t') in the x',t' system. 
The temporal dependence of the filter scale yields to the 
presence of a non commutation term which is also asso- 
ciated to the non stationary term. The trasformed non 
commutation terms, and the relevant approximations, re- 
leased by the unsteady and the convective terms, however 
cancel each other. Furthermore, the terms obtained from 
the divergence of the stress tensor in a system like (13T1 
21]) and the corresponding approximations ([53|) - ([5T|) are 
all Galileian invariants, which assures that the variable 
scale filtered equations and their approximations are also 
such. 

Another general implication, linked to the presence of 
a finite domain, is that the boundary conditions for the 
filtered variables should be different from those for the 
unfiltcred variables. The problem of wall boundary con- 
ditions for the filtered field could be treated with this pro- 
cedure by adopting one of the classical approximated con- 
ditions, which rely on the introduction of a special sub- 
grid model, that is, the wall model, which is apt to rep- 
resent the inner layer dynamics and which puts the first 
grid point inside the logarithmic laver[l5) ■ [l~6| ■ [j~7T| ' [l~8j . 
It could also be treated by placing grid points well in- 
side the viscous sublayer to resolve the near-wall dynam- 
ics and by assuming no slip and impermeability bound- 
ary conditions. It should be recalled that the latter 
conditions, which in theory should not be used for fil- 
tered velocities, introduce an error of 0(A 2 ), indepen- 
dently of the filter shape. Their use requires the sub- 
grid model, which should represent the non homogeneous 
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and anisotropic structure of the viscous and buffer lay- 
ers, to be altered in the inner region. The employement 
of anisotropic models based on a tensorial turbulent vis- 
cosity would be opportune, see Horiuti (1990) [19} . Carati 
and Cabot (1996) t20j, the review monography by Sagaut 
(2001, Chap. 5.3)[a] and also the differential angular mo- 
mentum model (Iovieno and Tordella, 2002) [2lj, which, 
being based on the representation of the turbulent vis- 
cosity through the moment of momentum vector, is well 
suited to assume an anisotropic formulation. 
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III. NUMERICAL TESTS 



A. A priori tests on the turbulent channel flow 

In this section a set of a priori tests is presented, which 
provides information on the field distribution of the non 
commutation terms, their relevant approximations and 
their ratios with respect to the physical terms from which 
they arise. The data correlating the approximated and 
exact non commutation terms have been determined by 
filtering the direct numerical simulation of the turbulent 
plane channel flow at Re T — 180, as performed by Al- 
fonsi et al. (1998) [lj and Passoni et al. (1999) and 
at Re T = 590, as performed by Moser et al. (1999) [Hj]. 
The longitudinal momentum balance, which implies a 
zero pressure non commutation term, is considered. Rep- 
etition on two different grid levels has been performed. 
All the data here presented have been averaged over an 
interval of 1.2 revolution times. 

Figure [T] shows the distributions of the exact non com- 
mutation terms and of their approximation according to 
the present procedure, at two filtering levels. The convec- 
tion term of the longitudinal momentum balance, its grid 
and subgrid scale decomposition and the diffusive term 
are shown in parts (a), (b) and (c), respectively. The 
filter <$(x) = (Ax,ip(y)Ay,Az), with constants Ax, Ay 
and Az, varies along the transversal non dimensional 
y direction according to <f(y), where ip(y) £ [0,1], and 
y £ [—1, 1] is a function of at least class C 2 . The vari- 
ation of ip(y) has been laterally arranged (in 20% of the 
channel width along the walls) as follows 



<p{y) = 



tanha(y + 1) tanha(l — y) 
tanh 2 a 



a = 4, 



(38) 



where a is the parameter that controls the gradient of 
the filter scale at the wall. The non commutation terms 
on the first and second derivatives have been determined 
for such an anisotropic structure of the filter through the 
use of relations (|A.16[) and (IA.20p . The data in Fig. 1 (a, 
c) show that the present procedure yields, on average, for 
the Reynolds stress 
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FIG. 1. (a): ( ) Average values of the exact non commu- 
tation convection term and (— — — ) correspondent average 
values predicted by the procedure, (b): Noncommutation ap- 
proximated terms for the resolved Reynolds stresses (— . — 

. — ) and for the subgrid scale stresses ( ). (c): ( ) 

Average values of the noncommutation diffusive term and (— 
) correspondent values predicted by the procedure. 



0.09, 
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\C |a 

and for the viscous stress 



2A 



0.18, 
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where A = 4.92 • 1(T 2 . 

These figures have been obtained by neglecting the 
data that belong to the first 5% near the wall, where the 
numerical uncertainty due to the spatial discretization is 
high - especially as regard the exact non commutation 
terms computations - and where these results deteriorate 
by nearly 25%. 

The results of the numerical test at Re T = 590 as com- 
pared with those at Re T = 180, see Fig. [T] (a), show a 
good invariance of the procedure accuracy with respect 
to variations of the flow control parameter. 

It should be remarked, that a procedure capable to 
predict at worst, using the rather large value of S = 
4.92 • 10~ 2 , the 90% of the value of the non commutation 
terms may be considered accurate. In fact, the analysis of 
numerical errors in LES of Turbulence with cutoff in the 
inertial range (errors due to spatial discretization: finite- 
differencing errors and aliasing errors) shows that the re- 
sulting errors are very large, of the same order and even 
larger than the magnitude of the subgrid term over most 
of the wavenumber interval, for finite difference schemes 
up to eighth-order accurate, irrespective of the grid reso- 
lution (cf. Ghosal S., 1996, pp. 201-202)|22|. In such a 
general situation, it may be considered it a success that 
the procedure is capable of predicting nine-tenth of the 
value of the non commutation terms. The relevant aver- 
age error cannot spoil the overall numerical reliability of 
the simulations since it is about one order of magnitude 
lower than the errors due to the spatial discretation. 

A warning is necessary regarding the numerical compu- 
tation of the exact values of the non commutation terms. 
Direct computation through the definition is not recom- 
mended, as, even with the implementation of a numerical 
differentiation of the sixth order of accuracy, it artificially 
amplifies the fluctuations that are naturally present in 
the data field. The exact values of the non commutation 
terms should be correctly evaluated by using the integral 
representation of the derivatives 8/ 85 or d/dSj, such as 
([3]) in Sec. Ql] or (|A.3|I in the Appendix. No such numeri- 
cal problems affect the computation of the approximated 
non commutation terms. 

Figure [2] provides information on the relative impor- 
tance of the exact non commutation terms with respect 
to the physical terms which causes them. In part (a) 
the average ratio \C y ( (uv) ) /d y ( (uv)) | has been plotted for 
A = 4.92 -10~ 2 , A = 9.84 • 10~ 2 and also with an increase 
of the wall value of the stretching factor a = d y ip — 8. It 
can be seen, that close to the point where the Reynolds 
stress reaches its maximum (y + w 30), and the diver- 
gence therefore takes in the average very small values, 
ratio values as high as 0.38 are reached with the coarser 
grid. The comparison of these results with the results 
of Fig. [ljb) indicates that the average exact value of the 
convection non commutation term is of the same order 
than the average divergence of the subgrid stresses. Fur- 



thermore, it is interesting to observe that the doubling 
of the stretching factor a increases the relevance of the 
noncommutation terms in this region nearly as much as 
the doubling of the grid coarsening does, see in Fig. HJa) 
the near wall region where y + < 14. Figure Uta) also 
shows a positive comparison of the field integral value of 
\C y ({uv))/dj((uuj))\ given by Fureby and TaborQ with 
the distribution of the same ratio that has been yielded 
by the database used here. [TTI| ~ [i~2| 

With respect to the average ratio \C yy ((u))/dy y ((u})\, 
Fig. [D^b) yields maxima local values of about 100% for 
the coarser resolution, and about 60% for the finer one, 
close to where the relevant non commutation terms reach 
their local maxima near to the wall. Leaving aside local 
maximum values detached from the wall and relevant to 
the coarser resolution, this ratio, in the central part of 
the field, settles to constant lower values close to 0.2±0.1. 

B. Commutation error on an analytical solution 

As seen in Sec. Ill Al there are four type of non commu- 
tation terms in the variable scale filtered incompressible 
Navier-Stokes equations (I5T1) . (1521 . which are all source 
terms. If one limits the analysis of their influence on the 
flow solution to the determination of the field distribu- 
tion of the values they take with regards to the values 
taken by the original terms of the equations, one would 
mainly find the foreseen result that the non commutation 
terms are not negligible where the gradient of the filter 
scale is high. This is however not sufficient to under- 
stand the way, localized rather than extended, in which 
the commutation error affects the flow solution. 

For this purpose, it has here been considered useful to 
study the behaviour of an extremely simple flow model, a 
sort of conceptual model, which has two characteristics: 
- just one type of commutation source term is present, - 
its exact filtered solution, that is, the variable scale fil- 
tered solution not affected by the commutation error, is 
known and thus could be used as the reference solution. 
The second characteristics can only be obtained by fil- 
tering the exact solution of the unfiltered equation of the 
motion. 

On the other hand, to prove the efficiency of any given 
procedure for the correction of the commutation error, it 
is also necessary to know the exact filtered solution of a 
test flow, which, in turn, requires the knowledge of the 
exact flow unfiltered solution. 

Such a reference state cannot be found in a turbu- 
lent configuration of flow, for which no exact solution is 
available. Reference is therefore made to a laminar flow, 
which has an exact solution. One should note that, in 
such a case, the filtered equation of the motion, when 
the filter length is a function of the point but the com- 
mutation error correction is not considered, is identical 
to the unfiltered equation. 

The steady laminar incompressible channel flow has 
been selected as the test flow, since it has only one com- 



8 



A 




s 


0.2 


V 












A 






0.1 


_v 



















A = 4.92 l(T 2 ,a = 4 








A = 9.84 l(T 2 ,a = 4 








A = 4.92 10 , a = o 






Fureby and Tabor, Ref.[3] 

/ 





-0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 -0.0 




E r ' ^r,max 



C" / <U>" 



0.08 

0.06 
0.04 
0.02 


i 

-0.02 
-0.04 



-0.06 



0.4 



0.5 



20 40 60 80 100 120 140 160 180 



ij- 0.6 









A = 4.92 10 2 


o =4 


it / 

SI : 
!i : 






A = 9.84 10 1 
A = 4.92 10 1 


o =4 
o = 8 


III ., 
■j 


<\ /'V' ■ 
A;- • i \ 


\ /a\ 







oo '— ' ' ' ' ' 1 

-1.0 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 -0.0 

y 



FIG. 2. Average of the absolute value of the ratio between 
the exact non commutation terms and the physical correspon- 
dent terms in the motion equation, (a) Convection, ( ) 

\C y ((uv))/dj((uuj))\, (b): diffusion. Adimensionalization by 
means of channel semi- width and u T . 

mutation term, the diffusion one, see Q36p. and its solu- 
tion is analytically known. 

The non dimensional momentum equation for the 
steady incompressible channel flow is written as 

d 2 y u = Red xP = -^ , Vy 6(0,1) (39) 

where G is the modulus of the dimensional longitudinal 
pressure gradient and the adimensionalization is based 
on the channel width d, the gradient G and the density p 
(the non dimensional pressure gradient d x p results equal 
to —1). The boundary conditions are 

u(0) = 0, u(l) = (40) 

The corresponding filtered equation is 

d%(u) s (y)+Re = -C"((u) s ), Vye(0,l) (41) 

where C"({u)s) is set to zero to determine the solu- 
tion which neglects the commutation error and where 



FIG. 3. ( ) Filter scale across the channel (S — tp(y)A, 

a — 8, see (|44[) . Ill B); (- - - -) local values of the com- 
mutation term C" , see (|42[). referring to the diffusion term; 
( - ) R — E r j Ermax, relative commutation er- 
ror for the solution of the filtered non corrected equation 
(C" — 0) normalized with respect to the field peak value, 

E r = [{(u) — {u} eX act) /{u} exact]- 




0.0 0.1 0.2 0.3 0.4 0.5 

y 

FIG. 4. Absolute E a = (8/ 'Re) ((it) - (it) exact) ( ) and rela- 
tive E r — [((it) — (u) exact) /(u) exact] { ) error distributions 

of the filtered velocity distributions: A - without the com- 
mutation correction, B - with the commutation correction. 

Curve ( ): error recovery (E a (A) — E a (B))/E a (A) 

with the distance from the wall. The exact filtered velocity 
reference distribution is (u) exact = (Re/2)[y(l — y) — ^S 2 (y)]. 



C"((u)s) is approximated by 

c"({u) s ) = -^[((uy s ) 2S -( u y s ] 

6>(y) 2 +5(y)6"(y) 
26^) KW5>25 - WaJ(42J 

to determine the solution which accounts for the commu- 
tation error with the present procedure. The numerical 
solution of equation (|4"Tj) is determined by solving the 
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corresponding unsteady filtered equation 



dt(u) S -^- e d 2 y (u) 5 = l+-_r(y„i. Vy (0. I ) (4:-!) 



through a fourth order Runge-Kutta time integration 
scheme - carried out until the steady state is reached 
- coupled to a fourth order finite-difference discretiza- 
tion of the domain. The double level of average has been 
computed using a third order Hermitian quadrature for- 
mula. 

The filter of the exact solution gives the velocity dis- 
Re 1 
tribution (u) exact = ~ v) ~ 3 5 feOL which con- 

stitutes the reference on which the constrast between 
the commutation corrected filtered solution and the non 
corrected filtered solution is based. The filter scale 
varies along the transversal direction according to §{y) = 
(p(y)A, where <p(y) £ [0, 1], with A = 0.1, is a function 
of at least class C . The variation of ip(y) 



tanh2aytanh2a(l — y) 



tanh 2 a 



(44) 



has been laterally arranged in 20% of the channel width 
along the wall, setting the parameter a = 4, see Fig. [3] 

By contrasting the filtered exact solution, the absolute 
and relative errors relevant to the corrected and non cor- 
rected solutions are compared in Fig. [4] The correspond- 
ing error recovery, with the wall distance, is also shown in 
Figs. [2] and [3] It can be seen that the present procedure 
greatly reduces the commutation errors: in the central 
part of the flow, an almost full recovery is obtained. 

For a comparison of the distributions of the local value 
of the non commutation term and of the relative commu- 
tation error on the solution, see again Fig. [3] It is im- 
portant to observe that, by neglecting the commutation 
correction, the field results to be affected by a system- 
atic error not only in the region where the filter length 
varies, but also in the region where it is constant. This 
behaviour is due to the accumulation of errors on the ve- 
locity variable and its derivative, which is due to the lack 
of the two diffusion addenda (see (|4*Tj) , I02J)) that should 
enter the momentum balance equation. Even through 
these terms are significantly different from zero in a lim- 
ited portion of the flow they affect the entire field to a 
great extent. In the central part of the flow, where the 
non commutation term is very small (Fig. [3]) , the relative 
commutation error results to be of the same order as the 
peak value of the field, while the absolute error reaches 
its maximum value (see again Fig.U). 



IV. CONCLUSIONS 

A procedure to explicitly insert the correction terms in 
order to counteract the commutation error associated to 
the use of a variable filter scale in the filtered equations of 
motion is here presented. With this procedure it is possi- 
ble to directly compensate for the commutation error on 



the filtered field. The procedure uses volume average fil- 
tering, but more general filter operators are also possible. 
Both isotropic and fully anisotropic filtering configura- 
tions are considered. Approximated commutation terms, 
with an accuracy of the fourth order in the filter width, 
are inserted into the motion equations, which do not in- 
crease their differential order. The difficulties related to 
the addition of further boundary conditions are therefore 
avoided. The proposed representation of the commutator 
operators is based on truncated expansions in the filter 
width of finite difference approximations, that make use 
of a multilevel average operation. This fact suggests the 
joint use of the present procedure with subgrid models 
which need an explicit filtering of the equations of mo- 
tion, such as dynamic and mixed models. 

A set of a priori tests, with a plane channel flow DNS 
(Re T — 180) as a test field, proves the good correla- 
tion that the present procedure yields between the ap- 
proximate and "exact" non commutation terms. It also 
provide information on the relative importance of these 
terms with respect to the original physical terms. The in- 
fluence of the field resolution on the general non commu- 
tation term is confirmed to be 0(A 2 ). Asymptotically, 
the accuracy of the present approximation is expected 
to be <3(A 4 ). At the resolution levels corresponding to 
Re T = 180, these tests show a reduction of the abso- 
lute errors, after halting the reference filter scale, that is 
nearly 0(A 3 ). 

The filtering is a mathematical operation, that is not 
specific of the Navier-Stokes equations and is indepen- 
dent of the solution typology. When it is varied, the filter 
causes one kind of non commutation term for each dif- 
ferential term present in the equations of the motion. A 
test, for which the analytically exact (unfiltered and as a 
consequence filtered) solution is available, has been con- 
sidered to overcome the limitation of an analysis, based 
on the a posteriori determination of the relative order 
of magnitude of the non commutation terms, with re- 
gard to the original physical terms of the motion equa- 
tions, and to analyze the effects of the commutation er- 
ror on a flow solution. The chosen test flow is the two 
dimensional incompressible laminar channel flow, whose 
dynamics consist of the balance between the constant 
longitudinal pressure force and the viscous diffusion. In 
this case, only one type of commutation source term - 
diffusion - is present, and this is only of relative impor- 
tance in the lateral part of the flow, according to the 
filter gradient dependence. The error, due to the lack of 
a non commutation term in the motion equation, is how- 
ever also transferred to the central part of the flow. The 
result is a biased filtered velocity distribution where the 
relative error in the central part of the flow, where the 
filter gradient is zero, is of the same order of magnitude 
as the local maximum error of the field, which is situated 
at a distance from the wall of about 15% of the channel 
width. It has been shown that the present procedure can 
reduce the commutation error by one order of magnitude 
in the central part of the field. 
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where 1 = (1,1,1), the transformation rjj — 6j£j (with 
det (drji/d^k) = SiS 2 S 3 ) has been introduced and thus 

In such a situation, by virtue of the fact that 



APPENDIX: NON COMMUTATION 
APPROXIMATION FOR ANISOTROPIC 
FILTERS 

When the geometry of the flow domain requires 
stretching each direction independently, which implies 
S(x) = (<5i(x),<52(x),<5 3 (x)), it is opportune to adopt a 
class of integration volumes of the kind 



and an average operation for the variable 



d 

dS k 

3 



1 f c Q l At 
V l JV-, OXk 



(A.3) 



if) 8 = tT / /( x + V)dri = ±- I f{x, + S&M, 
V s Jv s V 1 J Vl 



(A. 2) it is obtained 
J 



3=1 



V- 1 /* rrK^J 0/.,. 1 f df s^dS, 1 f df , f> ,5/, v % 



5/ 



a*,- i 



1 Jv- 



Vi 



5/ 



5/ 



1 Jv- 



dxi 



3=1 



where <5*^ is the Kronecker unit tensor. As in Sec. [TT] 
[see relation Q] for the isotropic configuration, it results 
that the anisotropic filter of the derivative is a differential 
operator that acts on the filtered field: 

3=1 J 

The anisotropic non commutation term C z ', which is de- 
fined as 

<$(</>*) = <|£>*-^</>* (A.5) 

can now be represented through (|A.4[) as the sum of 
the products of the filter space derivatives and the fil- 
ter derivatives of the filtered variable: 

«vw = -EiJ- </)5 (A ' 6) 

3=1 3 



Proceeding in strict analogy with what has been done 
in Sec. |TTJ the anisotropic non commutation term can be 
approximated through second order centered finite dif- 
ferences: 

mf) s ) = - E (</w je , - tnss^+oisj) 

Again, using a Taylor expansion of the integrating func- 
tion in (|A.2[) . where only even derivatives appear since 
the domain of integration (|A.1[) is symmetric with re- 
spect to all the integration variables, expressions are 
obtained for quantities such as (f)g and {f)s±Sjep * n 
terms of /, and the filter width - which is now defined 
as Sj — A(pj(x),Vj, where A is a reference value for |<$| 
and < <fij(n) < l,Vj, x: 
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-0(A 4 



(A.. 



becomes 



where 



^</>«) = -E^4(«/>*>2^-^) 

j=i j 



(A.14) 



and 



dipj 1 



E 



3^3 



-^E%^ 2 W0)^e,+O(A 4 ) 

3=1 3 
«/>tf(x))^ 3 . ej - 

yE^wS+°( a4 ). ( a - 9 ) 



3=1 



since from (TOI) (g)x = g + 0(A 2 ), Vg. 

The basic approximation for the anisotropic non com- 
mutation term C' t is derived through these expansions as 



(A.15) 

For flow fields where the domain grid needs to be 
stretched along only one direction, say y and whose typ- 
ical examples are two dimensional wall-bounded flows, 
the last representation yields very simple approxima- 
tion formulae. By adopting the widely used notation 
8(x) = (Ax, (p(y)Ay, Az) in such a case with constants 
A:r, Ay and Az, the anisotropic approximation for the 
first derivative non commutation term results to be 

C 'y((f)8) = -%2j$) {W8hy{y)Ay ~ if)s) 

(A.16) 

The anisotropic non commutation term of the second 
derivatives, being defined by 



i=i 



((if) 8) 8+SjBj ~ ((/)*)«-<y,-e J 

(A.10) 



dipj 1 
'-^ dxi 2ipj 



E 



{iif)d)s+5 j e j _ (if) 8)5-8^ 

(A.ll) 



c u((f)8) = (g^)8 - faS(f)6> 
can be obtained as in Sec. Ql] [see (f23|) -(f25 j) ] 

CU(f)8) = -t&^(f)8 



(A.17) 



j=i 



d8i 



^ Fit- v 



( 9 2 



while the accuracy of the anisotropic first derivative com- 
mutation error and its approximation can be verified to 
be 



3=1 
3 



dxi d8jdx. 



(f)s) 



dxi dxi 



0(A 4 ) 



3 = 1 



Cl((f} 5 ) = C' i ((f} s ) + 0(A i ). 



(A.12) 



(A.13) 



Before deriving the approximated form for (|A.18j) . while 
wishing to mantain its fourth order of accuracy, it is use- 
ful to observe that the terms on the r.h.s. which contain 
d 2 

— or (f)x, k ^ j, by (IA.8I) . are of the same order as 
adjddk 

the remainder term. Therefore, they do not enter the 
approximation, which is 

3 



By filtering only in the j direction, expansion (|A.8[) be- 
comes 



cm)s) 



j2 5 ^ 5 i 



pa 



\2 r 



3 = 1 



25] 



(WshSjBj if) s 



((f) S hs jej ) = /(x) + 2AV|(x)0(x) + 0(A 4 ), 
which implies 

((f)8hSie-(f)8 = ((f)8)8 + 8 i e-((f)8)8-5^+0( A % C«/>5' = 



E , 

J=l J 



(di(f)s)25 j e j - d i(f)8\ (A.19) 



Again, in analogy with what has been done for the first 
derivatives, Eq. (|A. 16|i . the following is obtained 



2<p*(y) 



((f)sh<p(y)Ay ~ (f>8 



and, though keeping the same order of accuracy, the 
equivalent representation for both relations (jA.101 IA.11[) 



(dy(f)8hv(y)Ay- d y(f)8\ (A.20) 
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